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Abstract 

We establish a self-consistent variational framework that allows us to study 
numerically the non-equilibrium evolution of non-perturbative inhomogeneous 
field configurations including quantum backreaction effects. After discussing 
the practical merits and disadvantages of different approaches we provide a 
closed set of local and renormalizable update equations that determine the 
dynamical evolution of inhomogeneous condensates and can be implemented 
numerically. These incorporate self-consistently the backreaction of quantum 
fluctuations and particle production. This program requires the solution of 
a self-consistent inhomogeneous problem to provide initial Cauchy data for 
the inhomogeneous condensates and Green's functions. We provide a simple 
solvable ansatz for such an initial value problem for the Sine Gordon and c/> 4 
quantum field theories in one spatial dimension. We compare exact known 
results of the Sine Gordon model to this simple ansatz. We also study the 
linear sigma model in the large N limit in three spatial dimensions as a mi- 
croscopic model for pion production in ultrarelativistic collisions. We provide 
a solvable self-consistent ansatz for the initial value problem with cylindrical 
symmetry. For this case we also obtain a closed set of local and renormalized 
update equations that can be numerically implemented. A novel phenomenon 
of spinodal instabilities and pion production arises as a result of a Klein para- 
dox for large amplitude inhomogeneous condensate configurations. 

11.10.-z,11.10.Wx,11.10.Ex,25.75.-q 
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I. INTRODUCTION 



The dynamical evolution of spatially inhomogeneous semiclassical field configurations 
along with quantum fluctuations that dress them play an important role in many relevant 
problems in high energy and condensed matter physics. Such situations are ubiquitous in the 
description of semiclassical processes such as baryogenesis via sphaleron decay [|],|] , bubble 
and droplet nucleation during supercooled phase transitions as it could be the case in QCD 
HID, the formation of topological defects [|]|| and the possibility of formation of disoriented 
chiral condensates [fZ| — 110|| - Localized high energy and large amplitude field configurations 
are of particular relevance for modelling the dynamics of relativistic heavy ion collisions as 
was proposed for example by ||TTJ] in terms of collisions of Yang-Mills wave-packets. In this 
case an approximation to the dynamics of high initial multiplicity states would be to provide 
an initial value problem with a particular initial condition on the field at some proper time 
surface. By proposing a large amplitude field configuration localized in space in the form of 
a semiclassical coherent state, the features associated with an initial state of large energy 
density and high multiplicity are captured. For example such field could be a strong color 
field or electric field whose relaxation will result in high multiplicity particle production fllT 



The dynamics of coherent structures and semiclassical field configurations is also of the- 
oretical and experimental relevance in condensed matter physics within the context of trans- 



port phenomena mediated by soliton excitations in quasi-one dimensional systems [12-14 



Although there have been several studies of real-time dynamics of classical field config- 



urations within the context of sphaleron decay M, Non-Abelian wave-packet collisions [11 



and coherent structures in scalar field theories JjlJ including classical fi~6Hl8| and semiclas- 
sical [19] studies of disoriented chiral condensates, there have been very few and limited 
attempts to incorporate consistently the effects of quantum fluctuations upon the dynam- 
ics of semiclassical inhomogeneous configurations and mainly in one spatial dimension [|20| . 
Whereas the dynamics of inhomogeneous field configurations had been studied in the small 
amplitude regime P2"| |, the important regime of large amplitude and non-perturbative field 
configurations is still not understood. 

Recently there has been tremendous progress in studying the non-equilibrium dynamics 
of high energy and high multiplicity quantum states including self-consistently quantum 
fluctuations in the case of spatial or spatial rapidity homogeneity to study particle production 
in strong fields, dynamics of phase transitions and relaxation and the dynamics of inflation 



21| - [ 35| . In particular these methods had been used to study the formation of disoriented 



chiral condensates, and pion production and relaxation during the chiral phase transition 
in low-energy phenomenological models [0,0 and more recently initial quantum states of 
large multiplicity distributions with localized momenta [P5fl . 

These studies rely on an implementation of variational and large N |[23|| - | 25[| approxi- 



mations generalized to strongly out of equilibrium situations. These approximations 



39fl are specially suited for studying non-perturbative phenomena and in particular the 



largeiVapproximation is amenable to consistent improvement |2"1 



Since these methods allow one to study the full quantum evolution in a non-perturbative 
manner they are specially tailored to studying the non-equilibrium evolution of initial inho- 
mogeneous quantum states of large energy including quantum back-reaction effects. 

This is specially important in the modelling of relativistic heavy ion collisions in a man- 
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ner that allows one to include the effects of quantum backreaction and non-perturbative 
particle production. A useful description of the dynamical evolution in relativistic heavy 
ion collisions is based on a hydrodynamic picture in terms of a fluid that is described in 
local thermodynamic equilibrium by a local energy momentum tensor and an equation of 
state p3| , fnf . Although such a description is not only physically compelling but also exper- 
imentally consistent at the present energies [|3"U42|1, it is at best phenomenological and not 
justified from a microscopic theory. The hydrodynamic description relies on the assumption 
that as a result of the collision a large amount of energy is deposited in a small region of 
space-time. Boost invariant hydrodynamics should emerge when the energy density 

is much larger than the typical hadronic masses, and is a result of hydrodynamic scaling 
when the resulting pancake width goes to zero with center of mass energy. In this scaling 
limit the relevant degrees of freedom for describing the expansion are the fluid proper time 
and fluid rapidity. 

Although trying to study the emergence of boost invariant hydrodynamics within the 
framework of QCD is at this stage an impossible task, the features of hydrodynamical 
evolution and particle production are rather robust and should be captured by a microscopic 
albeit semi-phenomenological theory. 

In this article we begin our study of the non-equilibrium dynamics of inhomogeneous 
field configurations of large amplitude including the back-reaction of quantum fluctuations 
self-consistently and non-perturbatively. The focus and goal of this program is to study the 
relaxation of a strongly out of equilibrium inhomogeneous initial state, the description of 
particle production and the emergence of a hydrodynamical description from an ab initio 
full quantum evolution in a microscopic, phenomenologically inspired quantum field theory. 

This article is devoted to setting up the variational self-consistent problem and to study- 
ing the subtle but important issues of renormalization. The implementation of variational 
or large N approximations to study the dynamics requires: i) to determine the Cauchy data 
for the inhomogeneous expectation value of the field as well as the relevant Green's functions 
that account for backreaction. Since both the variational and the large N approximations 
lead to a self-consistent condition, the Cauchy data must be determined self-consistently. 
The solution to this self-consistent initial value problem for inhomogeneous field configura- 
tions is typically an extremely difficult and numerically intensive task, ii) setting up a closed 
set of local renormalized update equations that can be numerically implemented in current 
computers. Both the locality as well as renormalizability aspects are extremely important: 
current computational capabilities cannot handle non-local update equations because they 
are memory intensive. The numerical evolution will be implemented on a space-time mesh 
and the results should be insensitive to the mesh size, this requires renormalized evolution 
equations, which is also important to improve numerical stability. This issue is clearly very 
important for the quantum backreaction because the Green's functions that are necessary 
have short distance singularities that must be renormalized. 

iii) Finally the numerical implementation of the update equations obtained in ii) starting 
from the initial data provided by the first step. 

In this work we seek to provide the first two steps to implement this program, leaving 
the numerical implementation to a forthcoming article. 

The coupled evolution equations that we find more suitable for a numerical approach 
in the self-consistent Hartree and large N approximations are summarized below in eqns. 
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We focus on studying the initial value problem and the update equations first in relevant 
cases in 1 + 1 dimensions not only because it is a simpler setting but also because it is a 
relevant problem not only in quantum field theory but also in condensed matter physics 
in the case of quasi-one dimensional systems. By treating these examples in one spatial 
dimension we determine a strategy to set up the self-consistent initial value problem that 
can be implemented for particular three dimensional geometries relevant to the description 
of 'colliding pancakes' in heavy ion collisions. This strategy is based on mapping the initial 
value problem (in terms of the variational parameters) onto a Schrodinger-like problem 
whose spectrum of eigenvectors and eigenvalues is exactly known. 

We then move on to the linear sigma model description of pion physics and implement 
a large N self-consistent approximation and the strategy developed in the one dimensional 
case to set up the Cauchy data. In this case we find some novel phenomena associated with 
spinodal instabilities for large amplitude field configurations. Finally we address the issue 
of renormalization and obtain the (semi) local update equations. 

The article is organized as follows: in section II we introduce the different formulations 
to implement the variational and Hartree approximations, discussing the advantages and 
disadvantages of each of the methods and determine the method that we find the most 
suitable for numerical implementation at this stage. 

Section III is devoted to a description of the initial value problem. Section IV presents a 
detailed analysis of the self-consistent initial value problem for the cases of the Sine-Gordon 
and 4 theories in 1 + 1 space-time dimensions. After analyzing the vacuum structure in 
the variational approximation we propose simple self-consistent Ansatzse for the Cauchy 
data for both topological (kink) and pulse-like (lumps) field configurations by establishing a 
correspondence between the initial value problem in terms of the variational parameters and 
the eigenvalue problem of a Schrodinger operator whose spectrum is exactly known. The 
renormalization aspects are studied in detail. These examples are relevant for studying the 
dynamics of solitons dressed by the quantum fluctuations in 1 + 1 dimensions and to study 
the relaxation of large amplitude field configurations via meson production. 

In section V we study the linear sigma model description of low energy pion phenomenol- 
ogy by implementing a large N approximation in 3 + 1 dimensions. We find that the leading 
order in the large N is similar to a Hartree factorization which allows us to use the formu- 
lation presented in section II. 

We use the results of the 1 + 1 dimensional case to set up self-consistent Cauchy data 
for an initial inhomogeneous configuration with cylindrical symmetry that could model the 
initial state after the collision of Lorentz contracted 'pancakes'. The subtle renormalization 
aspects are discussed in detail and we find a phenomenon akin to the Klein paradox that 
results in spinodal instabilities for large amplitude inhomogeneous configurations. 

In section VI the dynamical evolution equations are summarized in terms of fully renor- 
malized quantities and the initial data on all of the Green's functions that are needed for 
the update is determined. Section VII proposes a numerical strategy for an implementa- 
tion of the dynamics. We summarize our strategy for the present and future work in the 
conclusions. 
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II. NON-EQUILIBRIUM DYNAMICS: METHODS 



There are several alternative ways to establish a variational self-consistent formulation 
of the dynamical evolution, and although all of them are equivalent, some particular imple- 
mentation may present some practical advantages over others. 

Thus we begin the program of setting up the variational self-consistent formulation by 
first studying the different implementations, discussing their merits and drawbacks, and 
choosing a particular formulation that would lead to an economical numerical implementa- 
tion. 



A. Time Dependent Variational Principle 

The method for obtaining a time dependent variational approximation to the non- 
equilibrium dynamics in quantum field theory is to start with Dirac's variational principle 
for obtaining the functional Schrodinger equation. 

This variational approach is implemented by proposing Gaussian wave functionals 
[BUI, whose parameters are determined self-consistently. This approach leads to 



the time dependent Hartree approximation, and becomes exact in the case of a scalar field 



in the vector representation of the O(N) group in the large N limit f24l , |25 



Dirac's variational principle proposes to extremize the effective action 

d 

dt<y\i—-H\%>; (vm = 1 (2.1 

and leads directly to the Schrodinger equation: 

6T = 0^{i^--H}\*>=0 (2.2) 

In the Schrodinger (ip) representation, the wave-functionals are given by 

t] =< v?|^ > 
and the functional Schrodinger equation becomes 

i^*[M = J d d x^-\j^ + \v V {x)W^x) + V\^nvA (2.3) 

Where we have specified d-spatial dimensions. In later sections we will address d = 1 and 
d = 3 separately in specific examples. 

In this representation the variational approach is implemented by proposing a Gaussian 



or Hartree trial wave functional P4p6| ~ H39| : 

* v [<p,t] = N{t)exp{-J d d x J d d y{(<p(x)-<f) c (x,t))K:(xJ-,t)(rty)-My,t))} + 

i f d d x Tr c {x,t){ip{x) - (f> c {x,t))\ (2.4) 



y; t) = G^(x, y, t)/4 - iE(x, y, t) (2.5) 



5 



where in (|2.5|) we have written the kernel K,(x, y; t) explicitly in terms of its real G and 
imaginary E parts. The real quantities ip c ; G ; X ; 7r c are variational and have the following 
meaning: 

c (x,t) =< W v \<p(x)\9 v > ; n c (x,t) =< V v \-i6/5<p(x)\9 v > 
G(x,y,t) =< V v \(p(x)(p{y)\V v > -<p c {x,t) c (y,t) 

< % v \tp(x,t) \ . > -<f> c {x,t) n c (y,t) = 2 [ d d z G(x,z;t)E(z,y;t) + ^-5 3 (x-y) 

i duly, t) J 2 



(2.6) 



The sub index „ stands for variational. 

In terms of these parameters the effective action for the trial wave functional becomes 

T(0 C , tt c , G, £)= J dt < V v \id/dt - H\V V > 

= J dtd d x L e ( x ,t) ?^j^ + J dtd d xd d yZ(x,y-,t) dG{X gt y]t) ^ 

-Jdt<H> (2.7) 

where 

<H> = [ d d x{\ii 2 + 2[SGS](x,x) + l -G-\x,x) + \iy V f 
J 2 o 2 

+1 lim V x V v G(x, y)+ < V[<p(x)] >} . (2.8) 

< H > is a constant and is a first integral of the motion. The equations that follow by 
varying the effective action with respect to the variational parameters are: 

.. . _ 2 . , . d < V > 
7r c (x,t)= V <p c {x,t) [x,t); 

0(f>c 

(j) c (x,t)= 7T c (x,t) 

G(x, y; t)= 2 J d d z[Z(x, z)G(z, y) + G(x, z)E{z, y)] 
t(x,y,t)= J d d z[-2Z(x,z)i;(z,y) + ^G- 1 (x,z)G- 1 (z,y)} 

+[\^l-^^{m^-v) (2-9) 

These equations form a closed system and are the Time Dependent Hartree (Fock) equa- 
tions (TDHF). They have the drawback that they are non-local and involve not only the 
Green's function G(x, y\ t) but also its inverse, furthermore the renormalizability aspects of 
the dynamical equations are not transparent in this formulation. A numerical implemen- 
tation of an initial value problem must provide Cauchy data for the variational parameters 
and implement a numerical update procedure that yields finite quantities, i.e. is free of short 
distance divergences. 



6 



There is an equilibrium by product of this method in the form of the variational effective 
potential for the case of static homogeneous configurations which is given by p6,39 



V eff ((f) c ) = - < V\H\V > ; <tf|y>|tf>=0 c ; <tf|tf>=l (2.10) 

with Q being the quantization volume and < H > is given in equation (|2.8|) . This effective 
potential is useful to establish the phase structure of the particular theory in this Gaussian 
variational approximation. 



B. Method of Equal- time Green's Functions 

The update equations for the kernels that appear in the Time Dependent Hartree varia- 
tional wave function are non-local as displayed by the set of equations (|2.9| ). We can trade 
this non-locality by increasing the number of equations. This is achieved by considering the 
first order update equations for the equal time Green's functions in terms of the Heisenberg 
field operator (p and its canonically conjugate momentum tt = (p. Consider the evolution 
equation of the following quantities, 

(f) c (x, t) =< <p(x , t) > ; ir c (x,t) =< <p(x,t) > 
G(x,y,t) =< ip(x,t)<p(y,t) > -(j) c {x : t)(j) c (y,t) 
D(x,y,t) =< <p(x,t)n(y,t) > -cj> c (x,t)Tr c (y,t) 

D(x, y, t) =< tt(x, t)(p(y, t) > -(f> c (y, t)ir c (x, t) = D(y, x, t) - i5 d (x - y) 
K(x,y,t) =< (p(x,t)(p(y,t) > --K c (x,t)Tc c (y,t). (2.11) 

We can relate D D and K to G and S as follows: 

D(x, y; t) + D(y, x; t)= 4 J d d zG(x, z; t)E(z, y; t) 

D(x, y; t) + D(y, x, t)= 4 J d d zE(x, z; t)G(z, y; t) 

(2.12) 



K(x,y;t) = ^G-\x,y;t)+4: J d d z 1 d d z 2 E{x, z x ;t)G{z x , z 2 ;t)E(z 2 ,y;t) (2.13) 

This last equation, eq. ( |2.13| ) shows that K is not an independent quantity. 

All of the expectation values are with respect to the initial Gaussian trial wave-functional 

or density matrix. 

The time derivatives of operators are obtained from the Heisenberg equations with the 
Hamiltonian 

H = fd d x {\^xf + U/if ■ W V + + V[tp(x)]} (2.14) 

and the canonical commutation relations between the field (p and its canonical conjugate 
momentum n. 
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These lead to the Heisenberg operator equations of motion: 



ir(x, t)= (p(x, t) 

7r(x, t)= V 2 ip(x, t) - dV/dip(x, t) 



The evolution equation for </> c is as given before 



d 2 cj> c (x,t) + < 



dV 



>= 



d(f(x, t) 

and we obtain the following coupled equations for the kernels G, D, D and K: 

G{x, m t) = D(x, y; t) + D(x, y; t) 



D(x, y- 1) = K{x, y; t) + V 2 v G{x, y-t)- < {<p{x, t) - c (f, t)} 



dV 



d<p(y,t) 



> 



(2.15) 
(2.16) 

(2.17) 
(2.18) 



K(x, y; t)= V 2 x D(x, y; t) + V 2 y D(x , y; t) 

dV 

- < {<p(x,t) - 4> c {x,t)} 



d(p(y,t) 



> + < 



dV 



d(p(x, t) 



{<p(y,t) - My,t)} > 



+ < 



dV dV 

{<p(x, t) - (f) c (x, t)} > + < {ip(y, t) - <j> e (y, t)}- 



d<p(y,t) 



dcp(x, t) 



> 



(2.19) 



The set of eqns. ( p,16| - pT^ ) is equivalent to the TDHF equations (|2.9| ). 

If we evaluate all the equations of motion values using a Gaussian Density Matrix, the 
equations for <f) c , G, D and K close, which is the usual truncation of the Dyson equations 
in a mean field approximation. 

We would like to show that the equation for K guarantees in the update that the quantity 

N(x, y, t) = K(x, y; t) - \g~\x, y-t) - 4 J d d Zl d d z 2 Z(x, z x \ t)G( Zl , z 2 ; t)Z(z 2 , y- 1) (2.20) 



is conserved in time. In the spatially homogeneous case we showed 01 that N(x, y, t) was 
related to the initial number of particles present at t — 0. If we are in a pure state such 
as one described by a Gaussian wave function, < a^a >= at t = as will be discussed 
below. A nonzero value of iV can be obtained if we start with a Gaussian density matrix as 



discussed in [27 



First one can show explicitly, for the field theories discussed in this article, using the 
equations of motion for G and S that for K defined by eq.( |2.13| ) K obeys eq. ( |2.19| ) 
so that the equation for K is automatically satisfied if we solve the TDHF equations ( [2.9|) . 



Conversely, we can rewrite the equation for iV ( 2.20 ) in terms of D and D using the relation: 

E(x, V,t) = ~J dzG-\x, z; t){D(z, y; t) + D(y, z; t)} = - J d d z{D{x, z; t) + D(z, x; 

(2.21) 
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Using the equations of motion for G D D and K one can then show that =rr = 0. Thus if 
we choose our initial condition on K to satisfy the constraint: 

K(x,y;t = 0) = ^G- 1 (x,y;t = 0) 

+ - A J d d Zl d d z 2 {D(x, z x - 0) + D(zi, x- O)}^- 1 ^!, z 2 ; 0){D(z 2 , y; 0) + D(y, z 2 ; 0)} 

(2.22) 



then the conservation law (which is equivalent to eq. ( 2.1£ ) guarantees that this constraint 
is observed at all times. 

There is a definite advantage in these update equations over the TDHF ones, as we 
will see in specific examples below. The equal time Green's function equations lead to a 
spatially local update (apart from gradients), which is quite an advance over the previous 
(TDHF) case which required non local updates as well as matrix inversions and therefore 
is computationally more intensive. The small price one pays is to increase the number of 
variables by one. 



C. Covariant approach to mean field theory 

The covariant approach to mean field theory, arises from a self consistent truncation of 
the Schwinger Dyson hierarchy of equations for the N-point Green's functions to just the first 
two equations: namely an equation for the one and two particle connected Green's Functions. 
This can be accomplished directly in terms of the Dyson equations (as we do here), or by 
introducing a quadratic constraint field in a path integral approach and performing a loop 
expansion with respect to the constraint field ]23| . 



Starting from the Lagrangian 

L = J d d xdt {-df.ipd^if - V[<f}} (2.23) 

we want to obtain the equations of motion for the expectation values 

(j) c (x, t) = Trpcp(x, t) =< tp(x, t) > . (2.24) 
and the two time Wightman functions : G <,> (x,y;t,t') given by 

G > (x, y; t, t') = Tr[pip(x, t)(p(y, t')) - <f) c (x, t)(j> c (y, t') 

G<(x , y- 1, f) = Tr[ PV (y, t)<p(x, t')} - cf> c {y, t)(f> c (x, t') (2.25) 

Expectation values are in the initial density matrix p. The various Green's Functions can be 
constructed from the two Wightman functions. In the mean field approximation, the exact 
equation of motion for C is: 

9 2 c (x,t)+< 7 f^4->=O (2.26) 
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which is approximated by replacing the exact expectation value by the expectation value 
using the Gaussian trial wave-functional. Making the same approximation to the exact 
equation for the two point function one obtains for the Green's function: 

{d 2 + < ( ^^ )2 >}G(x,y;t,t') = 5 d (x - y)5 c (t - f) (2.27) 

Here because the boundary conditions pertain to an initial value problem, the delta function 
is defined on the Closed Time Path (CTP) contour of the Schwinger-Keldysh formulation 
44|. Namely 



fc(t _ = !«^Q (2 . 28) 

where 

The expectation value is obtained using the Gaussian density matrix. Since the expectation 
values obtained in the Gaussian trial state depend only on (fi c and G the system of coupled 
Green's function equations truncates at this level. The causal Green's function is related to 
the Wightman functions defined above as follows: 

G(x, y, t, t') = Q c (t, t')G>(x, y, t, t') + Q c (t', t)G<(x, y, t, t>). (2.29) 

In the covariant approach one needs two time information about the Green's function. Thus 
if we directly solve this equation, starting from the initial data at t — t' — 0, it would 
be computationally much more storage intensive than considering the set of coupled equal 
time Green's functions. One way to solve this Green's function equation and simultaneously 
handle the mass renormalization problem is to convert the partial differential equation into 
an integral equation. Defining M 2 (x, t) via 

MH X ,t) jfp% > 
{d<p(x,t)) 2 

and calling m 2 the mass in the Hartree approximation in the vacuum sector, then we have 
that 

G(x, y, t, t') = G (x, y, t, t') - J d d z J dt"G (x, z, t, t")AM 2 (z, t")G(z, y, t", t') (2.30) 
where 

AM 2 (x,t)) = M 2 (x,t)-m 2 r 

and 



dt = dt- dt 

Q,C+ J0,C- 



The CTP formalism [44] guarantees causality so that t" < {t,t'}. Solving this integral 
equation numerically is quite storage intensive. Here Go is the free field theory CTP Green's 
function for a particle with mass m 2 .. Integral equations of this type and their numerical 
solution are discussed in [[43|]. 
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To circumvent these storage problems, one can use another strategy for updating G 
which relies on finding a complete set of solutions for the quantum fluctuation field out of 
which the Green's function can be constructed. That is if we have a complete orthonormal 
set of solutions ipk(%,t) to the equation: 



{d 2 + M 2 {x,t)}M%,t) = 0, 



where the scalar product is defined by 
< f,g>=ij d d x 



r(M)^M- 9 (z, t ) s/ *^ f) 



Of 



Of 



(2.31) 



(2.32) 



then we can expand the fluctuation field tp (where <p — <fi c + -0) as follows: 

^ *) = ZX a * i^k(,x, t) + 4 ipl(x, t)} (2.33) 

k 

with the canonical momentum ir(x,t) given by: 

tt(x, t) = J2i a k d t ipk(x, t) + a], d t tpl(x, t)} (2.34) 
k 

(for continuum modes J2k = J (f^)- 

The mode functions satisfy the orthogonality conditions: 

< ip k , ipi >= 5(k -I) ; < ^l, ipi >= . 

Assuming the usual commutation relations for the creation and annihilation operators: 

[ajfe, aj] = 8ja\ [a k , a{\ = [a{, a}) = . (2.35) 

(where 5^ = (2n) d 5 d (k — I) for continuum modes) leads to the closure or completeness 
condition: 



Of 



Of 



8(x - y) . 



[2.36) 



One can verify using the equation of motion Q2.31 ) for the modes, that if the modes are 
orthogonal at t — that they stay orthogonal at all times. 

For the case of spatial homogeneity, it is quite simple to find a complete set of orthonormal 
solutions to the classical field problem at all times 
for spatially homogeneous problems 



{d 2 + M 2 (t)}Mx,t) =0 

Letting 

Mx,t)=fk(t) e tkx 
we find the fk obey the ordinary differential equation: 



The mode equations become 

(2.37) 
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f k + {k 2 + M 2 (t)}f k = 



In terms of the f k one then has that the fluctuation field can be constructed at all times as: 

$(x, t) = £{a fc /*(*) e tkx + 4 f* k (t) e- te } (2.38) 

k 

The closure relationship is then satisfied as long as f k obey the Wronskian condition at 
t = 0: 

f dfk r*dfk _ . 

Renormalization constraints require that for large values of k the modes f k differ in a pre- 
scribed fashion from the free field ones. The renormalization can be performed by studying 
the WKB expansion of the mode functions and isolating the divergences as discussed in 



several references: [^,^]. The initial values of fk and fk are usually chosen so that the 
initial state is either in an adiabatic vacuum state or a thermal state with respect to the 
adiabatic vacuum as discussed in |45].|29|. 



All the non-equilibrium Green's functions, can be constructed in terms of the Wightman 
function 

G > (x,y,t,t') =<<0(£,t)^(j/,O > (2.39) 

in the usual fashion. 

In the spatially inhomogeneous case, the mode functions obey the partial differential 
equation: 

{d 2 + M 2 (x,t)}ij k (x,t) = (2.40) 

M 2 (x, t) is determined from the same modes that solve this equation, thus it is necessary, in 
order to get started at t = to find a self consistent choice for M 2 (x,t = 0) and the mode 
functions. This will allow us to determine the initial conditions for the partial differential 
equations for each mode. This aspect of the problem will be addressed in the following 
sections. The second problem in this approach is to find efficient symplectic integration 
methods which will allow us to simultaneously update many partial differential equations, 
each of which is similar to updating the time dependent Schrodinger equation. The way 
that the initial conditions are handled is by assuming at t — that we can write: 

Mx,t) = , 2 , 1/2 g k (x); g* k = g- k 

Then the closure relationship becomes the completeness relation: 

Y,g* k (x)g k (y) = 6(x-y) (2.41) 

k 

and the g k obey a Schrodinger equation. 

{_V 2 + M 2 (x, t = 0)}g k (x) = uo 2 k g k (x). (2.42) 
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The covariant approach to mean field theory is often the first term in a systematic 
expansion of the full field theory, so that the advantage here is that one can determine the 
next order corrections and see how long the mean field approximation is valid. This has 
been carried out in a simple quantum mechanical example in ||43|| . The standard method 
for doing the systematic expansion is to start with the path integral and use closed time 



path boundary conditions to insure causality This is discussed in great detail in |24|,|25 



The advantage of having a mode function decomposition for the Green's functions is that 
subtractions necessary for renormalization can be made on a mode by mode basis, which 
dramatically improves the convergence of the mode sums. 



D. Operator and Closed Time Path Methods: 

A convenient implementation of the Gaussian variational approximation at the operator 
level begins with the Lagrangian formulation, shifting the field by its expectation value and 
implementing a Hartree-like factorization of non-linear interaction terms in order to make 
the Lagrangian linear plus quadratic. The linear terms in the fluctuation are required to 
vanish, leading to an equation of motion for the expectation value, the remaining Lagrangian 
is a quadratic form in terms of the fluctuations. Passing to a Hamiltonian formulation for the 
operator fluctuations the Heisenberg equations of motion for the operators is obtained, and 
the expectation values are now understood in the Heisenberg representation in the initial 
state. 

Specifically, for a Lagrangian density of the form 

C[<p) = ±d ll <pd»<p-V(<p) (2.43) 

(where we have omitted the linear coupling to sources), the Gaussian approximation is 
implemented in the following steps: i) shift the field (p and its expectation value 

0(x, t) = <p c (x, t) + $(x, t) (2.44) 
= (2.45) 

ii) After this shift the potential has in general even and odd powers of the fluctuation ijj, 

V[ip c + il>] = V e [<p e ,ip] +V°[<p c ,i/>] (2.46) 

V e [cp c , = V[<p e ] + V 2 (^PP) (2-47) 

V^,]=J d -^l) (2.48) 

The Hartree factorization replaces the even and odd parts by 

V^[ipcM -> V[ip c ] + ^ 2 M 2 (x,t) +C(x,t) 

M 2 (x,t) = (^|^> (2.49) 
V^ c ,^}=i;V(x,t) 

V(f,t) = (^M| $= ^) (2.50) 
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Where the expectation values will be obtained self-consistently. The field independent 
term C is obtained by requesting that the expectation value of the potential equals the 
expectation value of its Hartree factorized form. 

Upon integration by parts in the resulting action, the linear term in the fluctuation ip is 
requested to vanish to ensure the condition ( |2.45| ), leading to the equation of motion for the 
expectation value 



dt 2 



VVc + V(x,t) = 



(2.51) 



Passing to the Hamiltonian in terms of the operators ip and their canonical momenta 7r 



2 



d d x 



7T 2 + (V^) 2 + M 2 (x , t)^ 2 + C(x, t) 
the Heisenberg equations of motion for the operators become 

= 



^-V 2 + M 2 (x,t) 



The Hartree factorization requires the quantity 

G(x,x;t) = (^ 2 (x,t)) 



[2.52) 



(2.53) 



(2.54) 



in which < • • • > stands for expectation value in the initial Gaussian wave-functional or 
density matrix since in the Heisenberg picture the operators evolve in time but the states 
do not. Since V and M 2 are expectation values of the first and second derivative of the 
potential as given by eqs. ( |2.49| , |2.50[) this implies a functional relation of both V ; M 2 on 
G(x,x;t) leading to a self-consistent condition. 

An alternative description of the non-equilibrium dynamics is based on the Keldysh or 
closed time path method |44j] in which the main ingredient is the time evolution of an initial 
density matrix, which could describe either a pure or mixed state 



p(t) = U{t,t )p(t )U~ 1 {t,t ) 



(2.55) 



where U(t, t ) is the time evolution operator, which is completely determined by the Hamil- 
tonian of the system, with boundary condition U(t ,t ) = 1. 

Non-equilibrium expectation values and correlation functions in this time evolved density 
matrix can be handily obtained via a path integral representation of the generating functional 



of non-equilibrium correlation functions along a contour in complex time 



< Dlf + Vip - e iJdtd d x{c[ v ,+ ,j+]-C[^,j-]} 



(2.56) 



where the fields and sources ^ ; refer to the forward (+) and backward (-) time branches, 
and functional derivatives with respect to the sources generate all the non-equilibrium corre- 



lation functions. The reader is referred to the ample literature for more details [44,28,25,31 



The Gaussian approximation is implemented by performing the shift of fields given by 
eq. (|2.44j) with the same <£> c on both branches since it is the expectation value and with the 
factorization given by eqs. ( |2 . 49| , p~50|) for the fields on both branches. 
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In this factorization scheme and after an integration by parts in the generating functional, 
the effective lagrangian density becomes 



£{ Vc , V] = C[cp c ] + \ {d^d^ - ^ 2 M 2 } + i> {-d^tp, - V} (2.57) 

The effective Lagrangian density necessary for the calculation of non-equilibrium correlation 
functions in the CTP-Keldysh scheme is given by £ e //[<£> + , cp~] = C[ip c , ip + ] — C[<p c , ip~]. Thus 
we see that this factorization scheme leads to a Gaussian generating functional in which the 
correlation functions are obtained self-consistently. Examples are described in detail below. 
The condition ( |2.45| ) is ensured by requiring the vanishing of the term linear in if). This 



is the tadpole equation |3l[ that leads to the equation of motion for the expectation value 
given by ( |2.51| ). 



An alternative manner to confirm the Hartree factorization given by eqs. fl2.46 , 2~4"T , 2.48 



in this scheme is to write in the original Lagrangian density the even and odd parts of the 



potential after the shift ( 2.44|) as 



V e [<p c , V] = V[ip c ] + ^ 2 M 2 (x, t) + 5 e (2.58) 
V°[ip c ^} = ^V(x,t) + 5° (2.59) 

and treat 5 e '° as 'counterterms'. The Gaussian approximation is completely determined by 
the one and two point functions, and these counterterms are obtained by requiring that 5° 
gives a vanishing contribution to the one point function and 5 e gives a vanishing contribution 
to the two point function. This requirement leads to the identifications of V and M 2 given 
in eqs. ( ^49l^50|) . 



In the CTP formulation the important non-equilibrium Green's functions are given by 
the time-ordered Q ++ , anti-time ordered Q and Wightman Q + ~ , Q~ + functions given by 

g ++ (x } t; x', t') = Oj + (x } t)ij + (x', t')} = (Tip(x, t)i){x\ t')) (2.60) 

C}—{x, t; x', t') = (ip-(x, t)ip-(x', t')) = (f^(x, t)^{x', t')) (2.61) 

g + -(x } t-x\t r ) = -(^+(x, t)<T(a;', t')> = -{ip(x',t')iP(x,t)) (2.62) 

G~ + {x, t] x', t') = -{il)-{x, t)^ + (x', t')) = -(if)(x, t)^{x', t')) (2.63) 

With T ; T the time and anti-time ordering operators respectively. 

Since the Hartree approximation is local in space and time V and M 2 are local in space 
and time, the self-consistent dynamics only involves the spatial coincidence limit of the 
following Wightman function 

G{x,x';t) = ($(f,t)$(f,t)) = -g + -(x,t;x',t) = ~g~ + (x, t; x', t) (2.64) 

The time evolution of this Green's function leads to a closed Hamiltonian system of 
equations equivalent to the method of equal time Green's functions presented above, eqns. 
( P-16f|TT9| ) with the identification of the self consistent mass M and "potential" V: 



G(x , x 1 ; t) = D(x, x*, t) + D(x, x', t) (2.65) 
D(x, x', t) = (V>(x, t)i)(x , t)) ; D(x, x\ t) = (ip(x, t)^(f , t)) = D^, x, t) + i5 d (x - x') (2.66) 
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D(x, x', t) = K{x, x', t) + (yl - M 2 (x, *)) G(x, x', t) 

D{x, x, t) = K(x, x\ t) + (v 2 x , - M 2 (x , tj) G{x, x, t) (2.67) 



K(x, x', t) = (ip(x, t)ip(x' } t)) (2.68) 
K(x, x, t) = (yl - M 2 (x , t)) D(x, x', t) + (y 2 x , - M 2 (x', t)) D(x, x\ t) (2.69) 

This system of equations is Hamiltonian and can be updated in time by providing initial 
Cauchy data on the three functions G, D, K at an initial time as a function of two variables 
x, x'. Since ultimately we only need the concidence limit to solve the self-consistent dynam- 
ical problem not all of the range of x' is necessary to be stored numerically and for a given 
value of x only values of x' in the immediate neighborhood are necessary (see discussion in 
section VII). 

Of all the different available formulations that we presented above, we prefer the closed- 
time path method with the Hartree factorization leading to the evolution equations ([2.65- 



|2.69| ) as these present an economical and implementable system of update equations. How- 
ever the main purpose of enumerating and analyzing the different alternatives is to provide 
a detailed assessment of the advantages and shortcomings of each method. 

However before deciding on a numerical implementation of this system of equations two 
important steps are needed, the first is to provide self-consistent initial data for G,D,K 
and second is to address the renormalization issues to guarantee that the dynamics will not 
depend on ultraviolet (short distance) cutoffs. We will address these two issues in definite 
and relevant 1 + 1 and 3 + 1 dimensional examples below. 



III. INITIAL VALUE PROBLEM 

In all of the methods detailed in the previous section in which a closed Hamiltonian 
system of equations for the Green's functions determines completely the dynamics, the first 
step towards a dynamical implementation is to provide the initial value problem, i.e. the 
Cauchy data at the initial time for the Green's functions. 

The time evolution will be completely determined upon providing initial conditions on 
ip c and the Green's functions. The self-consistent approach requires to solve self-consistently 
the initial condition problem at t — 0. Once the solution to the self-consistent initial value 
problem is obtained the self-consistent time evolution is obtained from the closed time path 
set of equations ( |2.65| - |2.6^ ) obtained in the previous section. 



Because the set of equations are self-consistent, to provide this initial data is tantamount 
to solve the initial value problem for the Green's functions self-consistently for an initial 
inhomogeneous classical field configuration ip c (x,t = 0), where we have chosen t = as the 
initial time. 

Our strategy for setting up the initial value problem and the Cauchy data for the Green's 
function is more transparently described by focusing on the Heisenberg equations ( p.53|) . At 



the initial time we expand the quantum fluctuation field ij) and its canonical conjugate 
momentum it = ip in terms of a complete set of functions as 

$(2,t = 0) = J ^={a n ^ k (x)+ay*M} (3.1) 



16 



n($,t = 0)=i>frt = 0) = -i J d d k ^ {a s M^ ~ 4^} ( 3 ' 2 ) 

with ; at canonical annihilation and destruction operators. The mode functions ip k {x) 
are chosen to be the complete set of solutions of the following eigenvalue problem at the 
initial time 

'-u)\- V 2 + M 2 (x,t = 0)]il) k (x) = (3.3) 

In terms of these modes and the expansion of the quantum field and its canonical momentum 
given by eq. (|3.1|j3.2j ) the Hamiltonian (|2.52 ) at the initial time (t = 0) is given by 

H(t = 0) = J d d k u k (a\a k + 1/2) + / d d x C{x, t = 0) (3.4) 

Having found the complete set of mode functions, the initial values of G(x, x'; t = 
0) ; D(x,x';t = 0) ; K(x,x',t = 0) are completely determined by determining the ini- 
tial correlations {a\a k i) ; {dk^w) m the initial state. The most physically reasonable choice 
is given by 

(4 a k') = N k h,k' ; ( a k a k >) = (3.5) 

This choice corresponds to an initial state that is a density matrix which is diagonal in 
the basis of the number operator a\ak with occupation number N k . In particular for N k = 
the initial state is a pure state and corresponds to the Fock vacuum annihilated by the a k 
and is given by the instantaneous ground state of the Hamiltonian ( |2.52| ) at t — 0. 

At this stage we can make contact with the variational wave-functional by expanding the 
operator ip in the Schrodinger representation in terms of the complete set of states ^(x) at 
t = 

ip(x,t = 0) = J d d k r] k ip k (x) ; r)l = rj- k (3.6) 

where he choose the mode functions in such a way that ip k (x) = ip_ k (x) and the condition 
on the coefficients r\ k in ( p.6[ ) is a result of the hermiticity of the field ip(x). A pure state 
wavefunctional at the initial time can now be written as 

K k ,y = J d d xd d x> /C(f,f ) M#) 
n c ,k = / d d x n c (x) ^ k (x) (3.7) 



The choice leading to the initial correlations ( |3.5| ) corresponds to a wave-functional which 
is the ground state for all of the harmonic modes with frequencies u>k and 7r c = 0, i.e. 

K Kk , = ^ 5 k ^ ; ^k = (3.8) 

This choice leads to the identification of the real and imaginary parts of the kernel 

G k \i = 2uk5k,k' ; £fc,fc' = (3.9) 



For an analogous approach see |46 
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IV. 1 + 1 DIMENSIONS: KINKS AND LUMPS 



We begin by implementing the strategy to set up the initial Cauchy data as described 
in the previous section in two relevant examples in 1 + 1 dimensions. Not only does 1 + 1 
dimensions provide a simpler setting within which to gain experience and intuition, but also it 
is a relevant example for statistical mechanics and condensed matter physics since both Sine- 
Gordon and 4 theories are models for quasi-one dimensional charge density wave systems 
p~2|-|T4| . In these systems transport phenomena is mediated by coherent and collective 



excitations such as solitons, breathers and pulse-like excitations known as polarons |12Hl4 
Furthermore there is a large body of results on these quantum field theories |47] 
we can use as a yardstick against which to compare our approximations 



El that 



A. Sine-Gordon field theory 

The exactly integrable Sine-Gordon field theory provides a good example to test the reli- 
ability of variational methods since the quantum S-matrix is exactly known for all particles 
as well as the expectation values of local operators on asymptotic states. This theory is also 
known to have a phase transition as the coupling constant is increased [f|!|.|>2[ . 

The Sine-Gordon field theory is described by the Hamiltonian: 



H 



dx[—7i(x, ty 



i(V^,t)) 2 -pcos^] 



(4-1) 



Following the steps described in detail in the previous sections, we find that the Hartree 
factorization for this problem leads to 



V{x,t) 



M 2 (x,t) = a e-^ G{x > X]t) cos/fy c (x,t) 
G(x,x;t) = (ip 2 (x,t)) 



(4.2) 

(4.3) 
(4.4) 



The equations of motion that determine the dynamics in the Gaussian approximation 
thus become 



0Vc 


&<Pc 


dt 2 


dx 2 


~d 2 


d 2 


dt 2 


dx 2 



+ 



e 2 



1 sin j3ip c = 



-%-G(x,s;t) 



cos {3<p c 



tp(x,t) = 



(4.5) 



with the self-consistency condition given by eq. ( |4.4| ). 



1. 'Vacuum' sector: 



As mentioned above, the Sine-Gordon field theory has a phase transition at a critical 



value of the coupling constant |49|j52| . We will now show that this Gaussian approximation 
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reproduces this results of the Sine-Gordon field theory, giving us confidence in the approxi- 
mation. For this we consider the vacuum sector with a homogeneous value of the expectation 
value ip c = modulo 2tt. In this simple case the Heisenberg equations are solved in terms 
of plane waves, leading to 



G (x,y) = -!- J dkg (k)e 



ik(x—y) 



In the vacuum sector we find for the expectation value of the Hamiltonian density in the 
Gaussian state 

where we have introduced an upper momentum cutoff in the fc-integral. We now define the 
renormalized mass parameter 

«r = Mj^f'*" (4-8) 

with Gq = (tp 2 ) in the vacuum state and k 2 an arbitrary renormalization scale. The energy 
density once subtracted at the renormalization scale k 2 is given by 

,2 



Htf) - H(n 2 ) = ^ - k 2 ) - (4.9) 
Choosing the subtraction point to be 



in " (3 2 k 2 ' 



then the extremum of the energy is at 



k 2 = a r 



2 2 
jl = K . 



Since the second derivative of the energy is given by: 



<m K(^-i){^i p2/ ^- 2 (4.io) 



(d^ 2 ) 2 8nK 2y 8n /y K 2 
we find that a stable ground state exists only for 



f- < 1 (4-11) 
Sir 



which is in accord with the result of Coleman H49| . For (3 2 > 8tt the Sine-Gordon theory 
undergoes a Kosterlitz-Thouless phase transition [^]. Whereas the theory is superrenormal- 
izable for f3 2 < 87r it becomes renormalizable at (3 2 > 8tt since the cosine operator becomes 
marginal at this value of (3 [p^j . The region j3 2 > 8tt must be studied within a double 
expansion in a r and 5 = f3 2 — 8n [52]. 

We will restrict our study of the Sine-Gordon field theory to the superrenormalizable 
regime i.e. (3 2 < 8ir since the Gaussian approximation cannot describe correctly the physics 
for f3 2 > 8tt H|. 
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2. Initial value problem: Self-consistent Kinks and Lumps 



At the classical level, the Sine-Gordon model has moving kink solutions |47 



(x, t) — — tan 1 |exp 



Q 



(1-^2)1/ 



(4.12) 



where go is the position of the kink at time t — and v its velocity. The dependence of the 
solution on x,t reflects the translational and Lorentz invariance of the theory. The energy 



of this moving kink is pj| , [50 



E 



a 



1 



,2U/2 



(4.13) 



Thus the classical (c-number) kink is a relativistically invariant model of an extended body 
in one spatial dimension. 

For a kink solution at rest we have: 



cos (3(j)(x)= 1 



cosh aox 



sin (34>(x) 



2 sinh aox 
cosh 2 a x 



(4.14) 



We seek to set up an initial value problem that can include as self-consistent solutions 
kink-like configurations to be identified as a classical kink dressed by the quantum correc- 
tions. As described in a previous section, this is achieved by seeking self-consistent solutions 
of the following eigenvalue problem 



i 2 

dx 



2-tL-G(x,x) 



cos /3ip c (x) 







G(x, x) = G(x, x) — Gq 



(4.15) 
(4.16) 



where /i 2 ; Go are given by eq. (|4.6|) and with ip c a classical configuration with a kink-like 
profile to be found self-consistently. In the presence of a kink-like configuration the quantum 
fluctuations asymptotically far away from the kink will be similar to vacuum fluctuations 
since tp c {x) — * 2tt far away from the kink, therefore G(x, x) — Gq — > as x — » oo. This leads 
to the effective mass of the excitations to be given by fi. 

The spectrum of harmonic quantum fluctuations around soliton solutions of the Sine- 
Gordon theory is obtained from a Schrodinger like equation with a a potential of the form 



v{x) 



m 



n(n + 1) 
cosh 2 \mx] 



(4.17) 



where m is the (asymptotic) mass of the fluctuations and for the Sine-Gordon model n = 1. 
Potentials of this form are reflectionless and integrable and their spectrum is exactly known 
Therefore we propose a self-consistent potential for the mode functions of the form 



2 -81 
fj, e 2 



d( - x ' x) cos/3^ c (x) = M 2 — 



2m 2 



cosh \mx\ 



(4.18) 
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where the parameters M 2 ; m? are identified with the asymptotic mass of the quanta and 
the width of the kink or lump configuration and will be self-consistently determined from 
the parameters of the theory /i 2 ; (3. 

We want to stress that eq. (|4.18|) is an additional ansatz, independent of the variational 
approximation defined by eqs.( |4.4j )-( |4.5| ). We remark that although Eqs. ( [4.4| )-( p~5|) admit 
much more general solutions, the simple ansatz (|4.18|) is introduced to facilitate an analytic 
self-consistent solution to the initial value problem. Eqs.( |4.4| )-( |4.5| ) in general require a 
numerical solution (see sec. VII) which is computationally intensive. 

With this ansatz the eigenvalue problem becomes 



d 2 

dx~ 2+W ' 



m 2 + 



2m 2 



cosh \mx] 



i, k {x) = ; Wl 



m 



M 



(4.19) 



In terms of the solutions of this eigenvalue problem we will choose the initial Gaussian state 
to be that of the ground state for the Harmonic oscillators of frequency uo k . 

The solutions to this eigenvalue problem are available in the literature |47|.|53 



1 



2 coshlmxl 



W 2 



ikx 



[ik — mtanhfmx]] ; W k 



= 

k 2 + m 2 



M 2 -m 2 



k 2 + M 2 



(4.20) 
(4.21) 



We see that the mass of the asymptotic particle states around the classical configuration 
is given by M allowing us to identify at once this parameter with the (renormalized) mass 
parameter in the Sine Gordon equation 



M 



(4.22) 



Requiring the frequency of the bound state ( |4.20| ) to be positive restricts M 2 > 



m 



the continuum modes asymptotically are phase shifted plane waves. The Heisenberg field 
operators and its canonical conjugate momentum at the initial time are expanded as 



$(x,t = 0) 



(oq + oj) 1p (x) + 



1 



V2cJfc 



fc(x, t = 0) = -iJ -^(a - aj) ^o(x) - i ^ 



[ a k tp k (x) - a\ ipl(x) 



(4.23) 
(4.24) 



As mentioned above the initial state is chosen to be the ground state for the Harmonic 
oscillators of frequency u k , i.e. this state is the vacuum of all the annihilation operators. 
We find 



G(x,x) = G + 1 
cosh [mx\ 

Z7T 



m 



(4.25) 
(4.26) 



The ansatz ( |4.18| ) leads to the following consistency equations 
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M 2 = fj 2 

I3 2 F[ V ] 

cos (3(p c (x) = e 2cosh2 i m:E ) 



1 - 



2 V 2 



l + rf cosh 2 [ 



mx\ 



(4.27) 



The first equation is simply the statement that the asymptotic states far away from the 
soliton must describe quanta of mass fi, which has already been argued above in eq.( |4.22j ). 
For a topological kink solution we must impose that <p c (x — > — oo) — > ; ip c (x — > +oo) — > 2n. 
This requires that cos f3ip c (x — > ±oo) — > 1 but it must cross zero twice when the argument 
is 7i / 2 ; 37r/2 and that the minimum value of the cosine is —1. These two conditions lead to 



r] 2 > 1 

1 

2 



In 



(4.28) 
(4.29) 



The second condition ( |4.29|) determines m 
in the interval 



y - v 

m(/3) and for all values of (3 it is constrained 



— <m 2 < a 2 
2 ~ ~ P 



(4.30) 



this is a rather mild constraint on the width of the soliton and confirms that in the weak 
coupling limit (3 — > m — > \i ~ cto as is expected. Within the region ( [4.30| ) there are 
solutions for arbitrary (3. 

This method also allows us to construct self-consistent solutions of the initial value 
problem for the case of non-topological lumps. For all values of the parameters (still with 
ji 2 > m 2 ) such that 



rf-l 
r] 2 + l 



< 1 



(4.31) 



the solution provides a classical field configuration whose asymptotic value resides in the 
same vacuum state and is clearly an excited state that will relax into mesons under time 
evolution. 

In this variational approximation, physical quantities, such as the soliton profile, the 
energy, etc, are analytic functions of rj 2 . It is then interesting to connect rj 2 with the SG 



coupling constant (3. In the weak coupling limit {(3 
precisely, eqs. ( (4.26|) and ( |4.29| ) yield, 



0), we see that 77 — > 00. More 



7/ 



2/3 



1 + 



2_(§; 

3tt 1 4 



2/3 



+ o(/? 4/3 ) 



Therefore, the variational approximation yields for weak coupling a power series in rj~ 2 ~ 
/5 4//3 . This should be compared with the exact results which can be expressed in power series 
of j3 2 (besides essential singularities at (3 = 0). In other words, the variational approximation 
does not reproduce the analytic properties in the coupling constant since a cubic-root cut in 
(3 is suggested and the weak coupling theory and the variational approximation cannot be 
simply matched in a series in (3. 
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B. Comparison with exact results in the quantum sine-Gordon theory 



The expectation values of some interesting operators are exactly known in one-particle 
(and multi-particle) states p^J53| . Among them is cos/50 which is just the trace of the 
energy-momentum tensor, 



TO 



2«p 



[1 — cos (3<f>(x)] . 



Let us consider (quantum) one-soliton states \9 > with rapidity 9. That is a physical 
state with energy and momentum: 

E = M cosh 9 , p = M sinh 9 

where M is the physical quantum soliton mass. 

The expectation value of Tjf(0) in such states can be written as, 



f{9) e 1 \mo)\e 2 > 



Here 9 = 9 x -9 2 



One has the exact formula 55 



where 



and 



F{9) 



ft 2 



M sinh 9 
'If sinh f| 



(4.32) 



(4.33) 
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j(7T-0^ 
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\x9 




sinh [ttx] sinh 




cosh 





(4.34) 



[Our normalization convention for the soliton states differs on a factor M from ref. [55]]. 

The integral in eq. ( |4.34|) yields the so-called minimal form factor [|54j . It can be written 
as a double-infinite product of Gamma functions. In this way one explicitly sees that F(9) 
is a meromorphic function of 9. The closest pole of F(9) to the real #-axis is at 9 = ±2i£. 

For large momentum transfer, the expectation value ( |4.32| ) exhibits Regge behavior; 



F{9) 



e»i 



4 ^ M 



iir ) 
'Til 



where t = AM 2 cosh 2 | stands for the momentum transfer. 



The variational expectation value of cos[/30(x)] [eq. (|4.18|) 1 depends on the coordinate x 
indicating that the variational soliton state is not a momentum eigenstate but a superposition 
of momentum eigenstates, 



>= J dp f(p) \p = Msinh# > 
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where \p > stands for quantum one-soliton states and f(p) is some function normalized such 
that 

J \f(p)\ 2 dp=l. (4.35) 

Therefore, 

< *|T£(x)|¥ > = / dpdp' fipYfip') e«^*> < P \T^0)\P > 

J — oo 

= / dpdp' F(8(p,p')) e'G-f > a /(pO /(p)* (4.36) 

J — oo 

where F(9) is the form factor ( [O^ ) and ^ 2 cosh0 = a/M 2 + pPy/M' 2 + p 2 - pp' . 
Taking the Fourier transform of eq.( 4.36j ) yields, 

/ + 00 Arn . f 
— e «* < ¥\T£(x)\* > = jdq f(q)* f(q + p) F{9 M ) 



where 9 p>q = Argsinh ^ — Argsinh^j- (4.37) 



We read from eq. ( |4.18| ) that 



9 / 2m 2 
<*|T;(:r)|*> wr =- M 2 



/3 2 y cosh 2 [ma;] 
Therefore, 

-OO Arn 0/n 

e lxp < mTZ(x)\* >. P 



OO 



2vr MV /? 2 sinh 



7Tp 



2m 



(4.38) 



We can match the exact expressions and the variational approximation in the weak 
coupling limit which is also a non-relativistic limit p,p' << M, 9 « 1 since the solitons 
becomes very massive for small (3. Then, Q v q « jL and 



[ +00 p e*» < flWI* > P =° / + °° ^ < ¥|2£(a0|* >™/=° 

J-oo 27T ^ J-oo 27T Z7T 



where we used eqs. (|4.33|) , (|4.35|) and ( |4.37| ) and that M = 8m//3 2 for the soliton mass in the 
classical (3 — > limit. 

We do not see for large 9 Regge behaviour in F var (9). Eq. Q4.38D suggests that F var {9) 
decreases exponentially for large momentum transfer. 

In summary, the simple self-consistent ansatz (|4.18 ) does not reproduce the exact on-shell 
sine- Gordon except in the weak coupling regime. 

Thus this initial configuration will undergo non-trivial time evolution to relax via meson 
emission to a stationary state of lowest energy compatible with the variational ansatz. 

Since for weak coupling the variational initial configuration does not differ much from the 
exact solution, we expect slow time evolution at weak coupling with little meson production. 
Clearly a deeper insight into the further evolution of these variational initial configurations 
require a numerical study of eqs.([O D - (|4.5|) [see sec. VII]. 
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C. ip 4 Field Theory 



The v? 4 quantum field theory is the simplest theory that presents a symmetry breaking 
phase transition with relevant low-energy phenomenology in 3 + 1 dimensions. It is also a 
model to describe kink-like excitations and the collective behavior of quasi-one dimensional 
charge density wave systems within the Landau-Ginzburg framework. The kink configu- 
rations in the broken symmetry phase are interpreted as domain walls between coexisting 
degenerate phases. We will study the 3 + 1 dimensional version in the form of the linear 
sigma model in the next section. In this section we study the simpler 1 + 1 dimensional 
model not only as a preliminary step towards the 3 + 1 dimensional implementation but 
also because this is relevant in statistical mechanics and condensed matter to describe non- 
perturbatively the quantum dynamics of transport in quasi-one- dimensional charge density 
wave systems |T2| |14|. 

The Lagrangian density is given by 



C 



1 



-dr<pd»<p + ^ 2 - Y 



(4.39) 



Upon shifting the field by the expectation value (p c , the Hartree factorization advocated in 
the previous section leads to the following equations of motion 
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A* 2 +9^1 + 3G(x,x;t) 
ji 2 + 3g(fl + 3G(x, x; t) 



<f c (x,t) = 
rj)(x,t) = 



G(x,x;t) = (tp 2 (x,t)) 



(4.40) 

(4.41) 
(4.42) 



1. The vacuum sector and the phase transition 

Before proceeding to set up the initial value problem for inhomogeneous time dependent 
configurations we will address the properties of the vacuum sector in the variational approx- 
imation. For this we look at space-time independent configurations of ip c and obtain the 
effective potential V e ff(<p c ) = (7i) with H the Hamiltonian density as introduced in section 
II. For this purpose it proves convenient to introduce the composite operator 

X(x,t) = -n 2 + 3g£ 2 (x,t) (4.43) 

whose expectation value in the Gaussian variational state is given by 

x (x, t) = -/i 2 + 3gip 2 c (x, t) + 3gG(x, x; t) (4.44) 

which is recognized as the space-time dependent effective mass of the quantum fluctuations. 

For a space time constant configuration corresponding to the Fock vacuum of the anni- 
hilation operators, we find 

G(x,x) = [ — * (4.45) 
V ' ; J 2tt + x 



25 



and leads to the self-consistent gap equation 



X 



-fi 2 + 3g<p 2 c + 3g J 



dk 1 
2^ VP + X 



(4.46) 



Since the Heisenberg equations ( |4.40| , [4~4~1~1) must be finite, the logarithmic divergence in the 
integral is absorbed into a mass renormalization. This polynomial theory is superrenormal- 
izable in 1 + 1 dimensions and mass renormalization renders the theory finite. 
The effective potential is given by 



eff 



< V >-- 



l ..2j.2 , 9 
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(4.47) 



which can be written in terms of the two variables (p c ; % (the composite operator effective 
potential [|36[j ) in the form 



V.//(*x) = ^-f^-^ + A] + i/ 



dk 



2n 2u(k) 



(4.48) 



This form of the effective potential is useful because dV e ff((p c , x)/&X gives the gap equation 
(HI)- 

As a function of two variables we must minimize the potential with respect to both to 
find the vacuum state. The minimum of the potential is obtained from the condition 



X<\> - 2g<\? = 0. 



(4.49) 



The value of <fi at the minimum is the vacuum expectation value v and the expectation value 
of the composite operator evaluated at this minimum is given by 



X(j>=V 



2gv 2 . 



(4.50) 



The gap equation at the minimum is given by 



2gv' 



-\i + 3g 



v 2 + 



dk 



2tt 2uq (k) 



(4.51) 



where u (k) = (k 2 + m 2 ) 1 ^ 2 . We can use this value to subtract the gap equation for arbitrary 
(p c , and the renormalized gap equation then becomes: 



3<? 



X 



X = 2gv 2 + 3<?(0 2 - v 2 ) - \og[Aj 



An 



mi 



(4.52) 



Integrating the renormalized expression for |^ with respect to x gives us the renormalized 
effective potential which is: 



X i ±2 „,2\ 



2^ 



X-m r + xlog( — 
X 



(4.53) 
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and we have chosen the potential to be zero at the minimum for convenience. As a function 
of the coupling g this effective potential describes a (weak) first order phase transition as 
shown explicitly below. 

It is convenient to now rescale everything in terms of the renormalized mass m 2 = 2gv 2 . 
we let 

9 -> m 2 r g; , \ -> ™lx 



so that at the minimum we now have: 

x = i; 



1 

2^ 



The effective potential is now given by: 

f w - h - f u>' - ^: 

m A r 2 2g 2 4<j^ 



' (X-1) 2 + 7^(X-1-Xl0gx)- (4-54) 



12g v/v ' 8tt 
This can be most easily plotted by using the relationship: 

2g 3g 4tt 

and parametrically plotting V and as a function of x- To find the position of the phase 
transition, this takes place when the value of the potential at the minimum at <p = has 
the same value as at the minimum at = v which we chose to be zero. Solving for this 
numerically we find that 

gj m 2 r = .3896tt = 1.224. (4.55) 

To show this behavior we plot the effective potential at the three values g/m 2 . = l,g c ,2 in 
Fig. 1. 

For g < g c the potential has a double well structure with degenerate minima at ip c = ±v, 
as the value of the (dimensionless) coupling is increased the minima shifts abrubtly to = 
when g > g c , these results are similar to those found previously . 



2. Inhomogeneous Initial value problem: Kinks and Lumps 



The harmonic fluctuations around a stationary classical kink configuration are solutions 
of a Schrodinger-like eigenvalue problem with a potential given by eq. (|4.17|) for n = 2 whose 
solutions are also known exactly The classical kink configuration corresponds to 



1 



cosh \mx] 



(4.56) 



with v the tree level vacuum expectation value. Thus we seek a self-consistent solution of 
the eigenvalue problem in the Gaussian variational approximation of the form 
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ip 2 c (x) = v 2 



a 



(^(x)) = (^ 2 (oc)> 



.4 



cosh [rax] 1 cosh [mx] 

We propose a variational ansatz similar to that of the Sine Gordon problem above 

2m 2 



(4.57) 



/i + 3g (p c + G(x, x) 



M 



cosh \mx] 



(4.58) 



where the variational parameters M ; m will be determined self-consistently. The logarith- 
mic divergence in < ip 2 (x) > is independent of x and is absorbed in mass renormalization, 
which is the only renormalization required in this superrenormalizable theory. Thus define 
the renormalized mass parameter 



(4.59) 



Requiring that <p c (x) asymptotically for large \x\ becomes an exact solution of the space-time 
dependent problem leads to the constraint 



- fx 2 R v + gv 3 
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Then the asymptotic solution of the eigenvalue problem leads to the identification 



M 



2/4 



(4.60) 



(4.61) 



which is the usual relationship between the mass of the asymptotic quanta and the vacuum 
expectation value in the broken symmetry phase. 

Since the ansatz ( 4.58Q is similar to that in the Sine Gordon case studied previously, the 
results for G(x, x) can be obtained from the previous case and is given by eq.( ^25[) with 



dk 



27T2VFTM 2 



(4.62) 



Finally the expressions ( |4.25| , fi~62"|) with the identifications given by (|4 . 5 9|J4 . 6 0| , |4 . 6 1| ) and the 
ansatz for the expectation value (|4.56| ) lead to the self-consistent solution 



a 



2m 2 



+ F[ V ] 



(4.63) 



with F[tj] given in ( |4. 251 , ^26]) . The only possibility for a topological kink profile is given by 

a = v 2 =r- pI(x) = v 2 tanh 2 [mx] (4.64) 

Requiring that the initial configuration <p c (x) have a topological kink profile, the self- 
consistency condition (|4.63j ) leads to the following transcendental equation for r\ 



3-r] 2 



3(1 + iflFfo] 



(4.65) 
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where we have used eqs. ([4.64jJ4.6Q , [On| ) and the expression for r\ given in eqs. ( |4.25| , f4~26"|) , 
the right hand side of ( [4.65 ) is recognized as the dimensionless coupling. It is clear that 
there are solutions for any value of the dimensionless coupling with < rj < 3. 

We note that this variational form for <f c (x) is not a stationary solution of the equation 



of motion (|4.40|) but must be taken as an initial configuration that will eventually evolve in 
time. 

This method also allows us to construct non-topological lump solutions with a < v 2 . 
These solutions correspond to initial inhomogeneous configurations that will relax via emis- 
sion of mesons. 

Although we have renormalized the theory at the initial time by mass renormalization, 
it is straightforward to prove that the same mass renormalization can be performed at all 
times. The main point is that the Heisenberg equations for the operator ip can be written 
in the form 



dt 2 



+ M 2 + 6~M 2 (x,t) 



if>(x,t) = 



(4.66) 



and the mass counterterm SM 2 (x, t) is a localized function of space. The Green's function 
can be calculated in perturbation theory in SM 2 (x, t) in terms of Feynman diagrams with 
free field propagators of mass M and mass insertions of 5M 2 (x, t). The zeroth order has 
a logarithmic divergence which is absorbed in mass renormalization and all higher order 
insertions are finite as long as the theory remains superrenormalizable, which in the case of 
Sine-Gordon requires f3 2 < 8tt. This is equivalent to obtaining the divergences in an operator 
product expansion. 



D. Variational approach vs. Collective Coordinates: 



A perturbative approach to quantization around a soliton configuration begins by obtain- 
ing the eigenmodes of the Schrodinger operator associated with the harmonic fluctuations 
around the classical kink [|5(|,[)7]]- [ 63| . There is an eigenmode with zero frequency associated 
with translational invariance 0,0- [B3| that must be treated non-perturbatively. This 
mode is isolated from the mode expansion and treated as a collective coordinate |50|j57 l- 
|k| (although see for a different approach) and perturbation theory is carried out for 
the non-zero frequency modes, since these can be studied in a small amplitude expansion. 
Collective coordinate quantization is complicated by the fact that the functional basis is 
not cartesian |57| , |60| , |6T| with complicated Jacobians and a prescription for normal ordering 
|58| , |59|j62] , |63[| . Although high order calculations of S- matrix elements [^8i^9j|62l , |63l and more 



recently of non-equilibrium properties [Bl| had been carried out, they represent an enormous 
technical feat. 

In order to make contact with the collective coordinate program of quantization and the 
treatment of zero modes, we recognize that in the self-consistent solution for the fluctuations, 
the mode that would correspond to the zero frequency mode in the perturbative approach 



is precisely ipo{x) with frequency cuq = V M 2 — m 2 (see |4.20|) , that indeed becomes zero 
for m = M. However this value of m leads to infrared divergences in our treatment (the 
variable rj given by eq. ( |4.25| ) diverges) , which are a reflection of the infrared divergences that 
appear in the perturbative approach. We notice, however, that self-consistent solutions exist 
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for a large range of the ratio m/Mas discussed above and there are no infrared pathologies 
associated with any value of m for which solutions are available. 

The apparent contradiction between the variational approach and the collective coor- 
dinate method is resolved when it is realized that i) the profile (p c {x) obtained from the 
self-consistent solution to the mode equations is not a stationary solution of the equation 
of motion including the quantum backreaction effects, ii) the zero mode in the perturba- 
tive treatment is a consequence of translational invariance, the change in the classical field 
profile under an infinitesimal shift of the origin is the zero mode of the fluctuation operator 
]5Il|,3^p7|- In the variational approach translational invariance is also present since 

the ansatz for the self-consistent potentials (|4.18| , |4.58|) could be written in terms of x — Xq 
with xq an arbitrary point (the center of the soliton or lump) and the mode functions will 
be functions of x — xq. However, unlike in the perturbative calculation in which an in- 
finitesimal translation of the center of the classical configuration leads to a zero mode whose 
wave-function is oc dtp c (x)/dx, in the self-consistent approach, such an invariance requires 
an infinitesimal translation of tp c and the self-consistent potential which involves G(x,x). 
Obviously the function obtained by an infinitesimal shift in xq from (p c is not a zero mode 
of the fluctuation operator, precisely because of the self-consistent potential that provides 
the quantum backreaction. 

Thus our interpretation of the resolution of the apparent contradiction is that the quan- 
tum backreaction effects that dress the soliton or lump solution also remove the degen- 
eracy associated with the zero mode. Although translational invariance is present, the 
self-consistent ansatz determines a profile centered at a particular point in space, shifting 
the center is tantamount to shifting the spatial dependence of all of the Green's functions 
and expectation values under which the expectation value of the Hamiltonian in the trial 
state is invariant. 

Thus this self-consistent variational method based on a trial ansatz, provides a non- 
perturbative manner of circumventing the infrared singularities associated with translational 
zero modes in the perturbative approach. 



V. 3 + 1 DIMENSIONS: INHOMOGENEOUS CHIRAL CONDENSATES 

One of our main motivations for studying inhomogeneous non-equilibrium problems is 
to provide a consistent framework to study the relaxation of high energy inhomogeneous 
configuration in heavy ion collisions via expansion into the vacuum. The main goal is to 
provide a microscopic non-equilibrium description of the dynamics leading to a hydrody- 
namic description and the scaling solutions. Of particular relevance is the study of the 
relaxation of inhomogeneous chiral condensates since this could be potential soft probes of 
the chiral phase transition in the form of enhanced pion distributions with distinct isospin 
correlations. In order to address these issues we must i) provide a model quantum field 
theory that incorporates the relevant low energy physics, ii) we must propose an initial state 
that incorporates the relevant geometry, leading to longitudinal or spherical expansion, iii) 
provide a consistent framework to study the microscopic dynamics. 

For the first step we invoke the linear model in terms of a sigma field and an isospin 
triplet of pseudoscalar pions. In order to maintain the phenomenology of PCAC in the 
Gaussian approximation we invoke the large N limit in which the sigma field and the pions 
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are in the vector representation of the O(N) group. For the initial inhomogeneous field con- 
figuration one should ideally propose either an inhomogeneous cylindrically symmetric field 
configuration in the sigma component (or possibly with an inhomogeneous component along 
a definite isospin direction if one wants to study a disoriented configuration) reminiscent of 
two pancakes or a disk configuration with an extension in the transverse direction of several 



The study of such a configuration including quantum effects as advocated in this work is 
an extremely complicated task and beyond the computational capabilities at this moment. 
We will be less ambitious and seek to understand the dynamics of a simpler configuration 
with spatial variation along only one coordinate, chosen to be the z-axis and translationally 
invariant in the perpendicular directions. Even this simpler scenario represents an advance 
in the treatment of semiclassical configurations plus backreaction. Lastly, the time evolution 
is to be studied with the methods detailed in the previous section. As in the simpler 1 + 1 
dimensional problems we begin by establishing the self-consistent initial value problem and 
Cauchy data on the necessary Green's functions and inhomogeneous configurations. We 
note that rather than invoking boost invariance that would imply specifying Cauchy data 
on a proper time hyper-surface as is the case in a hydrodynamical treatment, we emphasize 
that boost invariant hydrodynamics must arise from the full microscopic dynamics, and we 
set up the Cauchy data on an initial constant time slice (t — 0). 

The Lagrangian density for the linear sigma model with an explicitly symmetry breaking 
term that gives a mass to the pions is given by 



The linear sigma model description of low energy pion phenomenology and PCAC corre- 
sponds to iV = 4 with parameters fit to reproduce the pion decay constant F n = 93Mev, the 
pion mass M w ps 140Mev and the sigma mass (the value of y/s at which the phase shift in 
the pseudoscalar channel in n — it scattering goes through n/2) M a pa 600Mev. 

In order to implement the Gaussian approximation which is appropriate in the large N 
limit we begin by performing the shift along the sigma direction 



To leading order the large N-limit is implemented by considering a Hartree-like factorization 
(neglecting 1/N terms) 



fm. 



C = ■ &*$ - \fi $ ■ $ + ±($ ■ $) 2 - ha 




(5.1) 
(5.2) 




(5.3) 




(5.4) 
(5.5) 
(5.6) 
(5.7) 




and assuming 0(N — 1) invariance by writing 
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N-1 



7f(x, t) = i>(x, t) (1, 1, - • • , 1) => <7T ■ 7?) = N{<IP 2 ) 



(5.9) 



Alternatively the large N expansion is systematically implemented by introducing an 
auxiliary field p4| , ^5| , this latter procedure allows a consistent implementation beyond the 
leading order. To leading order the two methods are equivalent ||31|| . This factorization can 
also be understood by implementing the tadpole method described in section II. 

Thus to leading order, the large N approximation is similar to the Hartree variational 
approach described in section II. The Hartree and large N equations of motion are similar but 
not identical, they differ by a factor of three between the classical and quantum contributions 



The tadpole equation leads to the evolution equation for the expectation value 

J^- V 2 -fil + ^<p 2 c (x,t) + ^ 2 (x,t)) <p c (x,t)=h (5.10) 

which is the mean-field equation obtained via the auxiliary field method f2"i| , j2"5| . After 
the large N-Hartree factorization and the implementation of the tadpole condition, the 
Hamiltonian for the field ip becomes quadratic and the Heisenberg equation for this field 
becomes 



dt 2 



- V 2 + M 2 (f, t) 



A 



tp(x,t) = 
A 



M 2 (x,t) = -//g + y -<pl(x,t) + -(^ 2 (x,t)) 



(5.11) 



We will propose a configuration ip c that only depends on one spatial coordinate and time, 
i.e. (p c (z, t). This simplification will allow us to make progress in setting up the initial value 
problem and to make use of the methods developed for the 1 + 1 dimensional case. Thus for 
the initial value problem we will solve the eigenvalue problem ( |3.3| ) with the ansatz 



^(z) = F 2 + 



a 



cosh 2 [mz] 

M 2 (x) = -t4 + ^ 2 c (x, t=0) + ^(^ 2 (x, t = 0)) 



Ml 



2m 2 



cosh \mz] 



(5.12) 

(5.13) 
(5.14) 



Asymptotically the solutions should be phase shifted plane waves describing asymptotic 
pions quanta of mass M n and therefore asymptotically for x >> 1/m the PCAC relation 
M 2 F n = h is fulfilled with 



(5.15) 



is the renormalized pion mass. Translational invariance in the perpendicular directions allows 
us to write the mode functions solutions of the eigenvalue problem as 
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m 2 + 



2m 2 



il>k(x) 
dz 2 

W 2 k =ool-kl-Ml + m 2 



cosh [777.2:] 



Mz) = 



(5.16) 
(5.17) 
(5.18) 



We will take as initial state a Gaussian wavefunctional corresponding to the ground state 
of the harmonic oscillators with frequencies uj k - 

In terms of the frequencies W k the problem is now exactly the same as the 1 + 1 di- 
mensional problems that were solved self-consistently in the previous section leading to the 
eigenfunctions and frequencies given by 



Ak-x 



2 cosh [m z] 



k 2 , + M 2 - m 2 



[ik z — mtanh[m,2]] ; uo k = \J k 2 + M 2 



(5.19) 
(5.20) 



where again stability requires that M 2 > m 2 . In the perturbative quantization approach, the 
transverse modes associated with ipo(x) are the 'capillary' waves associated with massless 
fluctuations transverse to the interface with dispersion relation = \k\ as a consequence of 
the translational zero mode in the z-direction. 

At the initial time the Heisenberg field operator is expanded in terms of this complete 
set of mode functions. With a continuum normalization 

/d 3 k 1 f d 2 k 1 

ww^m { aCk)MS) + h ^ + J t^tts + h ^ (5 - 21) 
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{a (k)Mx) ~ h.c.} (5.22) 



^a(k)ipk{x) — h.c. | — % 

Taking the initial state to be the vacuum for all the annihilation operators we finally find 

(^(f)) = (^ 2 (oo))+/ 
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(5.23) 
(5.24) 



The first term is the usual coincidence limit of the two point function and has a quadratic 
and logarithmic divergence which as usual is absorbed in a mass renormalization. Introduc- 
ing an ultraviolet cutoff the integrals in I yield to 
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(5.25) 
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Finally the ansatz ( |5.12|J5.13| ) leads to the following self-consistent condition for the effective 
mass term in the mode equation 



2m 2 



Am 2 



cosh 2 [mz] 2 cosh 2 [mz] 



a + I 



(5.26) 



The logarithmic divergence is absorbed in coupling constant renormalization 

A 



1 



1 1 

A + 167T 2 



In 



(5.27) 



where we chose to renormalize at the scale M n . After renormalization, the self-consistent 
solution for the inhomogeneous expectation value becomes 



Mtt 



rn 



m 



F 2 cosh [mz] 



(5.28) 



This self-consistent solution for the initial value problem is not a stationary solution of 
the evolution equation for the expectation value. This field profile will evolve in time and 
dissipate into meson excitations. The dynamics of this process of relaxation is studied within 
the time dependent scheme described above and addressed in detail below. 

The mass renormalization and coupling constant renormalization make finite the effective 
space dependent mass term M 2 (x,t = 0), we now show that this renormalization renders 
the effective mass term finite self-consistently at all times. The proof is an implementation 
of the operator product expansion and goes as follows. Since the inhomogeneity is localized 
in space and asymptotically the expectation value ip c (x,t) reaches its vacuum value, we can 
write 



M 2 (x,t) = M 2 + SM 2 (x,t) 



(5.29) 



The renormalization is obtained by requiring that the equations of motion are finite which 
results in the requirement that the effective mass term M 2 (x, t) be finite for all values of the 
space time arguments. Since M 2 is the renormalized pion mass, renormalizability implies 
that 



SM 2 {x,t) 



A r 



G(x, x, t) — G + <fl(x, t) - (fl(oo) 



A 



R 



Gr(x, x, t) + (p 2 c (x, t) - (p 2 c {oo) 



(5.30) 



is finite. 

The divergences in G(x,x,t) can be obtained by implementing the operator product 
expansion or alternatively by treating the term SM 2 (x, t) as a localized perturbation, i.e. as 
a mass insertion in the Green's function. 

It is a straightforward exercise to expand the Green's function G(x, x, t) which is the 
coincidence limit of the Feynman propagator in a perturbative Feynman series. This is a 
tadpole diagram (a two-point function in the coincidence limit) in the presence of a back- 
ground mass insertion, the zeroth order (one loop) is simply Go = (ip 2 (oo)) given by eq. 
( |5.23| ) and which has been absorbed in the mass renormalization and is cancelled by the 
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regular (space-time independent) mass counterterm Go in ( |5.30D . The first order in SM to 
one loop is simply a mass insertion in the tadpole that has exactly the divergence of the 
one-loop contribution to coupling constant renormalization (actually the renormalization of 
the composite operator < ip 2 > which is equivalent to coupling renormalization in the large 
N limit). Since the divergence is independent of the momentum transfer into the loop and 
higher order insertions are finite by power counting we find the result 



G(x, x, t) — Gq 



In 



A 



5M 2 (x, t) + G R (x,x,t) 



(5.31) 



with Gr(x, x, t) being finite. This expression leads at once to the renormalization of coupling 
and mass just as at t = 0. Writing 



A = Z\\ 



16tt 2 



In 



A 



(5.32) 
(5.33) 



then the quantity 



& t) = 



G(x, x, t) —G + <fl(x, t) 



(5.34) 



is finite. A more elegant and systematic treatment of the renormalization along the lines 



described above can be offered following the analysis presented in reference [65]. The linear 
sigma model description of low energy pion phenomenology dictates that in order to repro- 
duce the values of F n w 93Mev, « 140Mev and the sigma mass ~ 600Mev, the value 
of the quartic coupling is extremely large Xr ~ 36. For such a large value of the coupling 
there is a Landau pole at a value of the cutoff A < 2 — 3 Gev requiring that this theory be 
understood as a low energy theory for scales well below this cutoff. 

In a weakly coupled theory the form of the initial self-consistent configuration ( |5.28| ) 
would allow to prepare initial inhomogeneous states for which ip c probes the false vacuum 
(f c ~ on length scales much larger than the Compton wavelength of the pion. This would 
be achieved by specifying a value m « M n such that 4m 2 / XrF 2 « 1. However in the 
strongly coupled linear sigma model and in the stable region for values of m < M n there 
is no range of the variational parameter m that would allow to probe the 'false vacuum' 
region <^ c ^ 0. The main reason for this constraint is that we restricted our study to the 
stable region for which M 2 > m 2 to avoid the instabilities associated with the bound state 
solutions given by eq. (|5.19|) . 



Spinodal Instabilities from Klein's paradox 



The nature of the instabilities for m 2 > M 2 is clear from the eigenvalue eq. ( |5-17| ) when 
written in the following form 



+ k 2 + Ml 



2m 2 



coshfmz] 



ip k (z) = uJ 2 k ip k (z) 



(5.35) 



35 



For m? > M 2 the potential penetrates the negative energy continuum in a region in space 
of width Az m 1/m and there will be a band of perpendicular momenta k±_ for which there 
will be instabilities arising from particle production. This situation is similar to that of the 
'Klein paradox' for free quantum fields in the presence of external potentials. In a region 
near the origin the potential is negative and strong enough to cross the threshold 2M n and 
produces particles. This instability towards particle production for the case m? > M 2 is 
also encoded in the typical square root in the expressions ( |5.25| ). 

For this case there is a band of wavectors k± for which there will be strong production 
of pions with soft perpendicular momentum. This mechanism of particle production is very 
similar to that of spinodal instabilities in the case of phase transitions, since the modes of 
the quantum field feel an inverted harmonic oscillator potential, i.e. an effective negative 
mass squared. We can deal with this interesting situation in the same manner as in previous 
treatment of the dynamics as an initial value problem in phase transitions with spinodal 
instabilities. 

The method consists of preparing the initial state wave functional for the unstable modes 
that is the ground state wave functional for upright harmonic oscillators, with reference 

frequencies u% = \Jk\ + /i 2 , with /i 2 an arbitrary positive mass parameter. This is the 
quantum field theory analog of preparing an initial Gaussian state, which is the ground 
state of a particular harmonic oscillator potential and evolving it in an inverted harmonic 
oscillator potential, which is the physical situation for spinodal instabilities. 

Therefore we consider the situation with m 2 > M 2 by expanding the field in terms of 



the eigensolutions ipk of G5.17Q as in Q5.21]) but with the frequencies 



«2 



A' 



k\ + /i 2 for k\ < A 2 

k\ - A 2 for k 2 ± > A 2 
m 2 - Ml > 



(5.36) 

(5.37) 
(5.38) 



and the initial state is annihilated by all of the annihilation operators. In this case we 
find the result 



<V> 2 (f,t = 0)) = (^ 2 (oo)> + / 



Air 

F y]=-zz 

87T 



In 



A 



+ F 



i + y 2 



m ' m 

-y 



4vr 2 



1 + -VT 



x 2 In 



i - vT 



i + \fi 



(5.39) 
(5.40) 

(5.41) 



with (^ 2 (oo)) given by QOBQ . 

The square roots and the logarithmic branch singularities are the hallmark of thresholds 
to particle production and provide an intuitive explanation for the spinodal instabilities for 
large amplitudes: the field profile probes the negative energy continuum in a localized region 
in space, resulting in particle production. 

Renormalization proceeds as in the previous case, and leads to the final form of the 
self-consistent field profile for the expectation value in the case in which m 2 > M 2 

,2 



m' m 



m 



F 2 cosh [mz] 



(5.42) 
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In this case by varying fi ; m we can provide an initial self-consistent inhomogenous 
state that probes field configurations with large amplitudes away from the vacuum (f c = F n 
in either direction. 

The band of unstable perpendicular wave- vectors will result in profuse production of soft 
pions through the spinodal instabilities with a likely enhancement of pions in the low pt 
region. Although these spinodal instabilities could also be present in the one dimensional 
problems studied in the previous section, in that case there is only one mode that is unstable. 
In addition, integrability forbids particle production in 1 + 1 dimensional sine-Gordon field 
theory. Physically the instability of one mode is not very interesting in the one-dimensional 
models where we are interested in the dynamics of soliton states. However it is certainly 
interesting and physically motivated in the sigma model in three spatial dimensions where we 
are interested in the dynamics of pion production from a coherent initial state. Furthermore, 
these spinodal instabilities associated with the Klein paradox phenomenon are not particular 
to the trial profiles studied here but will necessarily be a robust feature of inhomogeneous 
states for which large amplitudes of the order parameter can probe the negative energy 
continuum within a localized region in space, resulting in particle production. 

Fig. 2 shows the profile for (p c (z)/F n for a set of parameters that allows to probe near the 
false vacuum near the origin, and fig. 3 shows (p c (z)/F 7T for a different set of parameters that 
probes a region with large amplitudes away from the vacuum. Both configurations probe 
large amplitude field configurations away from the vacuum in a localized region in space of 
width Az pa 1/m. 



VI. DYNAMICAL EVOLUTION: SUMMARY OF RENORMALIZED UPDATE 
EQUATIONS AND INITIAL CONDITIONS 

Having provided the initial state which furnishes the Cauchy data for the proper Green's 
functions and addressed the renormalization issues we are now in position to provide the 
update algorithm to study the dynamical evolution. We begin by gathering all of the equa- 
tions of motion and initial conditions for ip c (x,t) and the Green's functions which are given 
byeqs. (p3T|j233|j235H2T69| ) 



ip c (x , t) = 7T c (x, t) (6.1) 
7r c (x,t) = VVc-V(^,G) (6.2) 
G{x, x; t, t) = D(x, x 1 , t) + D(x, x\ t) (6.3) 

D(x, x', t) + D(x, x', t) = 2K{x, x', t) + (v* + V 2 X , - M 2 (x, t) - M 2 (x', t)) G(x, x', t) (6.4) 
K(x, x' , t) = (V 2 X - M 2 (x, tj) D{x, x', t) + (V^ - M 2 (x' , tj) D(x, x' , t) (6.5) 

where we have made explicit that V depends on ip c (x,t) and G(x,x,t) and is given by the 
appropriate expressions in the 1 + 1 dimensional cases (see eqns. ( |4.2| , |4.3| , |4.40|J4.41|) ) or by 
M 2 (x,t) = M 2 (ip c ,G) in the 3 + 1 dimensional case in the large N limit (see eqn. ( |5.11| )). 
The initial conditions chosen for the dynamics of the expectation value are given by 



<p c (x, t = 0) = f C)SC (x) ; 7r c (x, t = 0) = 



(6.6) 
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where <p c> sc(x) is the self-consistent profile given by eqs. ( |4.27|) and ( }4.29|) for a topological 
profile or (}4.31|) for a lump profile in the Sine-Gordon case, or the profile given in eq. ([4.57 ) 
with ( |4.63| ) for the 4 theory in 1 + 1 dimensions. In the 3 + 1 dimensional sigma model in 
the large N, the profile for the expectation value is given by eq. (|5.28|) in the stable case, 
i.e. without spinodal instabilities and by eq. ( |5.42| ) in the spinodally unstable case. 

In order to update the equations we need to obtain the value of the kernels D, D, K at 
the initial time t = 0. This is achieved by using the Heisenberg equations of motion ( |2.53| ), 
the field expansion given by ( |3.1| , |3.2| ) in terms of the eigenmodes ip k of the differential eq. 
( fj.3| ) and the initial correlations Q3.5p. We must distinguish between the different cases 
without and with spinodal instabilities. 



A. Without spinodal instabilities 

In this case we find the following expressions for the kernels at the initial time: 

G{z, z', t = 0) = t^m. + f dk Mzm*) - ml + 1 dimensions (6 7) 

2coo J 27r 2uj k 

Since the eigenmodes are rather simple in form the integrals above can be carried out 
analytically in all cases, thus providing the initial Cauchy data for the Green's function. 
The initial data for the other kernels is determined by 

D(S, x', t = 0) = -- £ MWuffl = --/(x - x 1 ) (6.9) 

D(x,x',t = 0) = l -J2M?)r k (x') = l -5\x-x') (6.10) 
1 k 1 

K(x,x*,t = 0) = £ ^ k (x)r k (x') (6.11) 
k 1 

[V 2 X - M 2 (x, t = 0)] G(x, x', t = 0) = [V 2 X , - M 2 (x', t = 0)] G(x, x', t = 0) = 

-E^w:^) ( 6 - 12 ) 

k 

We see that whereas the individual kernels evaluated at the initial time are divergent, 
the combination of the kernels that enter in the update equations above are finite and lead 
to the following result for the initial value problem 

G(x, f, t = 0) = (6.13) 

D(x,x',t = 0) = D(x,x',t = 0) = (6.14) 
K(x,x,t = 0)=0 (6.15) 

If the expectation value (p c {x) were a solution of the time independent evolution equation 
for the expectation value, we would have found an equilibrium self-consistent solution of 
the quantum back-reaction problem. However, as we pointed out before, the self-consistent 
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profile for <p c (x) is not a static solution of the evolution equation and in the second step of 
the update, if c (x) will begin to evolve and with it the effective mass squared M 2 (x) for the 
quantum fluctuations. 



B. With spinodal instabilities in 3 + 1 dimensions: 



In this case case the Heisenberg field expansion is as given in eq.( p.21[ ) but with the 



frequencies u;° replaced by the frequencies given by eq. (|5.36H5.38|) . In this case we find the 
same expression for either the 1 + 1 or 3 + 1 dimensional cases where the frequencies in 
the first term in eq. (|6.8| ) are given by eqs.( |5.36| - p.38| ). 



G{x,x,t = 0)=0 (6.16) 

D(x,x',t = 0) + D(x,x',t = 0) = j — ^ 0(A 2 - k\) *-M*)1>q@) (6-17) 

K(x , f, t = 0) = . (6.18) 



Obviously the nonvanishing and positive expression ( |6.17|) evaluated at x — x' which 



corresponds to the second derivative of G with respect to time is a manifestation of the early 
exponential growth of the quantum fluctuations associated with the spinodal instabilities 
and particle production. 

The systems of update equations shown above displays explicitly that although the in- 
dividual kernels K, D have short distance singularities, these cancel in the combinations 
that are necessary for the update equations. Insofar as the Green's function G(x,x',t) is 
evaluated at separate spatial points it remains finite throughout the evolution since the 
initial value is finite and the update equations are in terms of finite quantities. However, 
once the coincidence limit is taken, quadratic and logarithmic divergences will appear in 
the short distance operator product expansion. From the time evolved value of G(x,x',t) 
we must now take the coincidence limit to construct the updated effective mass M 2 (x,t) 
and the short distance divergences must be renormalized. However our previous analysis 
of the short distance divergences determines unambiguously the manner to obtain the fully 
renormalized equations. In 1 + 1 dimensions the equations of motion are rendered finite by 
mass renormalization and the replacement 

r A dk 1 

G(z, z) -= G sub (z, z) = G(z, z) - / ——=== (6.19) 



2tt VP + M 2 

where A is the upper momentum cutoff. In three dimensions we can take the spatial coin- 
cidence limit of the Green's function in the finite quantity M 2 + 5M 2 (x,t) with SM 2 (x,t) 
given by the expression ( |5.34| ) . 



VII. NUMERICAL STRATEGY: 

Although implementing a full numerical study of the evolution equations with the self- 
consistent initial data is beyond the scope of this paper, we briefly describe how we will 
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implement numerically the evolution in a forthcoming study. The implementation is obvi- 
ously much simpler in the one dimensional case: one begins by discretizing space introducing 
a lattice cutoff a and the upper momentum cutoff becomes A = n/a, so that the integral for 
G can now be done at once leading to function of Ma. Since eventually we are interested in 
the coincidence limit of G(z, z') we do not need to update the functions of two variables z, z' 
for the whole range of the two variables, the non-locality of the update equations are mild 
only requiring the derivatives at separate points. Thus for fixed z we can simply update the 
equations for z' just a few lattice spacings away in such a way that the discretized derivative 
does not introduce coincidences. The coincidence limit will correspond to writing z' — ■> z + a 
and the renormalized Green's function is obtained by subtracting the value of G , in the 
limit when Ma « 1 the renormalized Green's function is insensitive to the short distance 
cutoff. 

In 3 + 1 dimensions the strategy consists of the discretization of space in a cubic lattice 
of equal spacing a along all the directions. Now the second derivatives needed in the update 
equations require storing the field configurations on 6 sites around the particular points 
x, x', again since ultimately we are interested in local field configurations and the spatial 
coincidence limit we do not need to store the configurations for all of the range of x' but only 
in the immediate neighborhood of the point x. For a lattice spacing a there corresponds an 
upper momemtum k max = (7i/a,7i/a,ir/a) and the momentum integrals needed for G are 
done with spherical symmetry and a radial upper momentum cutoff A = y/Sn / 'a which is 
the radius of the spherical Brillouin zone. The coincidence limit must be understood now as 
G(x, x + a) and for Ma « I the renormalized quantity M 2 (x, t) is insensitive to the short 
distance cutoff. 

As discussed in the previous section within the context of renormalization, the linear 
sigma model is a low energy theory with a cutoff A < 2 — 3 Gev and this restricts the value 
of the lattice spacing for a numerical analysis. There is, however, a much more stringent 
constraint on the sigma model description of low energy physics: this model does not include 
other important hadronic degrees of freedom such as the pseudoscalars rj ; rj' with masses 
~ 500 — 600 Mev and vector mesons (which determine pion electromagnetic form factors) 
with masses 700 — 800 Mev. These constraints suggest that the actual physical cutoff 
for the sigma model is not much bigger than A m 500 — 700 Mev resulting in a lattice 
spacing such that M n a « 0.2. Thus the combination of Landau pole arguments and low 
energy constraints on the dynamics require that a physically reasonable lattice spacing for 
a numerical treatment be of the order of ~ 0.2 fm, which will not result in strong short 
distance effects in the update of the dynamical equations. 

VIII. CONCLUSIONS: 

In this article we focused on establishing a consistent framework to study the non- 
equilibrium evolution of inhomogeneous non-perturbative field configurations including 
quantum backreaction effects. The main goal is to provide numerically implementable 
schemes to study the time evolution of large amplitude field configurations including the 
quantum backreaction associated with particle production on the relaxational dynamics. 

We proposed a self-consistent variational approach to study the non-perturbative aspects 
of the evolution of large amplitude inhomogeneous configuration, and analyzed the avail- 
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able frameworks to implement them. Since our main goal is to provide a practical method 
amenable to numerical implementation we highlighted in detail the advantages and short- 
comings of the different approaches. We recognized the main ingredients that are required 
for a practical numerical implementation of this program: i) a set of local dynamical update 
equations both in space and time, ii) a self-consistent solution of the initial value problem 
that provides the Cauchy data at the initial time, iii) renormalizability implemented at the 
level of the update equations, so as to ensure insensitivity to the lattice introduced in the 
numerical implementation. After analyzing several different manners to implement a vari- 
ational, self-consistent treatment of the non-equilibrium dynamics for these inhomogeneous 
problems, we identified a formulation that leads to a closed set of local update equations. 
We then studied in detail the renormalizability of this set of equations in both one and 
three spatial dimensions. The study of the dynamics of non-perturbative inhomogeneous 
field configurations in one spatial dimension is warranted by the interest in understanding 
the dynamics of soliton configurations in the presence of quantum excitations, but also as 
model systems to test these methods in simpler settings. In one spatial dimension we provide 
simple self-consistent ansatz for the initial Cauchy data in terms of exactly solvable poten- 
tials for the Sine-Gordon and c/> 4 model. In the Sine Gordon case the exact mass spectrum, 
S-matrix, form factors are known and this allowed us to compare this simple ansatz to exact 
quantum results. Numerical resolution of the full set of variational equations should provide 
better approximations. 

In three dimensions we focused on the linear sigma model in the large N limit as a 
model for relaxation of non-perturbative coherent pion configurations such as disoriented 
chiral condensates. In this case we provided a self-consistent solution for the initial Cauchy 
data in the case of cylindrical symmetry by extending the results from the one dimensional 
examples. 

We found that for large amplitude configurations a mechanism akin to the Klein paradox 
resulting in particle production is responsible for novel spinodal instabilities that will lead to 
profuse pion production at low p t through the relaxation of the inhomogeneous condensate. 

Furthermore we provided the full set of renormalized (semi)local update equations that 
can be implemented numerically and discussed in detail the strategy for a practical imple- 
mentation, which is the next step of the program and currently in progress. 
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FIG. 3. <p c (z)/Fn vs. mz for M w /m = 0.0019, n = , X R = 36 
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